Class 2
Zanim zaczniemy
0. Przydatne zasoby pomocowe w pracy z R - część
druga
- Podręcznik
tidyverse dplyrcheat-sheet- Tidy data cheat-sheet
- Automatyzacja
z
purrr - Ściąga
z obsługi dat i czasu w
lubridate - Podsumowania
danych z
overvieR - Jak
uczyć
R? - Ściąga z podstawowych pojęć statystycznych
- Wybór testu statystycznego
- Ściąga z wyrażeń regularnych
A. tidyverse oraz dplyr - “porządne
dane”
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.7 ✔ dplyr 1.0.9
## ✔ tidyr 1.2.0 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
Ekosystem tidyverse to niejako “państwo w państwie” -
rodzaj dialektu R w mocny sposób ingerujący w sposób, w
jaki obsługiwane są pewne stałe elementy języka. Podstawowa rola
tidyverse to stworzenie wspólnego i spójnego języka, w
którym można opisywać standardowe struktury danych (tibble,
dplyr), automatyzować funkcje i działania
(purrr) oraz standardowo wizualizować dane w oparciu o
gramatykę grafiki (ggplot2).
1. Filtrowanie danych: filter()
Filtrowanie odbywa się za pomocą par
zmienna OPERATOR warunek - podobnie jak przy filtrowaniu za
pomocą indeksów. dplyr upraszcza wszystko pozwalając nam
bezpośrednio odwoływać się do nazw zmiennych.
Przefiltrujmy dane tak, by wyświetlić tylko potomstwo
matki nr. R187528:
dane <- read.table(file = "BTdata.csv",
sep = ",", header = T,
stringsAsFactors = T)
filter(dane, dam == "R187528")## tarsus back dam fosternest hatchdate sex weight habitat bill_length
## 1 14.46 549.6392 R187528 A2602 45 Fem 9.7 forest 11.565
## 2 14.61 549.6629 R187528 A2602 45 Fem 10.0 forest 11.737
## 3 14.16 549.3603 R187528 A2602 45 Fem 9.4 park 11.378
## 4 14.08 552.7623 R187528 A2602 45 Fem 9.4 forest 11.261
## 5 14.08 552.5236 R187528 A2602 45 Fem 9.5 park 11.296
## 6 14.77 550.0366 R187528 A1302 45 Male 10.1 park 11.825
## 7 14.01 551.5308 R187528 A1302 45 Fem 9.7 park 11.225
## 8 15.07 549.1769 R187528 A1302 45 Male 10.1 park 12.100
## bill_depth
## 1 0.631
## 2 0.668
## 3 0.583
## 4 0.603
## 5 0.593
## 6 0.652
## 7 0.631
## 8 0.684
Czy wszystkie młode wychowywały się w jednym gnieździe
(fosternest)?
unique(filter(dane, dam == "R187528")$fosternest)## [1] A2602 A1302
## 104 Levels: A1002 A102 A1202 A1302 A1602 A1802 A18B02 A2202 A22B02 ... H702
Wybierzmy spośród tych piskląt tylko takie, których skok przekracza 14.50 mm, lub masa przekracza, lub jest równa, 9.5 g:
filter(dane, dam == "R187528" & (tarsus > 14.50 | weight >= 9.5))## tarsus back dam fosternest hatchdate sex weight habitat bill_length
## 1 14.46 549.6392 R187528 A2602 45 Fem 9.7 forest 11.565
## 2 14.61 549.6629 R187528 A2602 45 Fem 10.0 forest 11.737
## 3 14.08 552.5236 R187528 A2602 45 Fem 9.5 park 11.296
## 4 14.77 550.0366 R187528 A1302 45 Male 10.1 park 11.825
## 5 14.01 551.5308 R187528 A1302 45 Fem 9.7 park 11.225
## 6 15.07 549.1769 R187528 A1302 45 Male 10.1 park 12.100
## bill_depth
## 1 0.631
## 2 0.668
## 3 0.593
## 4 0.652
## 5 0.631
## 6 0.684
Zadanie Wybierz wszystkie obserwacje, które są
samicami i mają kolor płaszcza (back) o dominancie barwnej
550 nm (w zaokrągleniu w dół do liczb całkowitych; użyj
floor()), oraz wykluły się (hatchdate) przed
dniem 45 (15 maja). Ile jest takich obserwacji?
Spodziewany wynik
filter(dane, sex == "Fem" & floor(back) == 550 & hatchdate < 45)## tarsus back dam fosternest hatchdate sex weight habitat bill_length
## 1 14.31 550.0810 R187517 C602 44 Fem 9.5 park 11.434
## 2 15.22 550.1262 R187517 C602 44 Fem 10.1 park 12.206
## 3 14.61 550.2081 R187517 C602 44 Fem 9.5 forest 11.698
## 4 14.54 550.1580 R187517 C602 44 Fem 9.8 park 11.620
## 5 14.99 550.6216 R187517 C2202 44 Fem 10.3 park 12.042
## 6 14.69 550.4395 R187398 C602 44 Fem 10.1 forest 11.801
## 7 14.99 550.3093 R187517 C2202 44 Fem 10.4 park 12.021
## bill_depth
## 1 0.629
## 2 0.696
## 3 0.644
## 4 0.573
## 5 0.679
## 6 0.638
## 7 0.677
nrow(filter(dane, sex == "Fem" & floor(back) == 550 & hatchdate < 45))## [1] 7
Zadanie Zmodyfikuj powyższy filtr, by wyciągnąć samice o danej dominancie barwnej LUB te wyklute przed 15 maja. Ile jest takich obserwacji?
Spodziewany wynik
nrow(filter(dane, sex == "Fem" & (floor(back) == 550 | hatchdate < 45)))## [1] 90
\(\blacksquare\)
2. Sortowanie za pomocą arrange()
Posortujmy dane mające dominantę barwną 550 nm (jak wyżej) oraz wyklute przed 15 maja wg długości skoku.
dane2 <- filter(dane, floor(back) == 550 & hatchdate < 45)
arrange(dane2, tarsus)## tarsus back dam fosternest hatchdate sex weight habitat bill_length
## 1 14.16 550.4797 R187398 C602 44 Male 9.6 forest 11.295
## 2 14.27 550.5134 R046108 C402 44 UNK 9.6 forest 11.403
## 3 14.31 550.0810 R187517 C602 44 Fem 9.5 park 11.434
## 4 14.54 550.1580 R187517 C602 44 Fem 9.8 park 11.620
## 5 14.54 550.2464 R187517 C2202 44 Male 9.8 forest 11.644
## 6 14.61 550.2081 R187517 C602 44 Fem 9.5 forest 11.698
## 7 14.61 550.2466 R187517 C2202 44 Male 9.8 park 11.623
## 8 14.69 550.4395 R187398 C602 44 Fem 10.1 forest 11.801
## 9 14.69 550.3631 R187398 C602 44 Male 9.8 forest 11.719
## 10 14.99 550.6216 R187517 C2202 44 Fem 10.3 park 12.042
## 11 14.99 550.3093 R187517 C2202 44 Fem 10.4 park 12.021
## 12 15.22 550.1262 R187517 C602 44 Fem 10.1 park 12.206
## bill_depth
## 1 0.644
## 2 0.598
## 3 0.629
## 4 0.573
## 5 0.654
## 6 0.644
## 7 0.622
## 8 0.638
## 9 0.679
## 10 0.679
## 11 0.677
## 12 0.696
Zmodyfikujmy sortowanie tak, by uporządkowanie według długości skoku
odbywało się w obrębie każdej matki genetycznej
(dam).
arrange(dane2, dam, tarsus)## tarsus back dam fosternest hatchdate sex weight habitat bill_length
## 1 14.27 550.5134 R046108 C402 44 UNK 9.6 forest 11.403
## 2 14.16 550.4797 R187398 C602 44 Male 9.6 forest 11.295
## 3 14.69 550.4395 R187398 C602 44 Fem 10.1 forest 11.801
## 4 14.69 550.3631 R187398 C602 44 Male 9.8 forest 11.719
## 5 14.31 550.0810 R187517 C602 44 Fem 9.5 park 11.434
## 6 14.54 550.1580 R187517 C602 44 Fem 9.8 park 11.620
## 7 14.54 550.2464 R187517 C2202 44 Male 9.8 forest 11.644
## 8 14.61 550.2081 R187517 C602 44 Fem 9.5 forest 11.698
## 9 14.61 550.2466 R187517 C2202 44 Male 9.8 park 11.623
## 10 14.99 550.6216 R187517 C2202 44 Fem 10.3 park 12.042
## 11 14.99 550.3093 R187517 C2202 44 Fem 10.4 park 12.021
## 12 15.22 550.1262 R187517 C602 44 Fem 10.1 park 12.206
## bill_depth
## 1 0.598
## 2 0.644
## 3 0.638
## 4 0.679
## 5 0.629
## 6 0.573
## 7 0.654
## 8 0.644
## 9 0.622
## 10 0.679
## 11 0.677
## 12 0.696
Zadanie Posortuj cały zestaw danych według wagi ciała, a następnie według gniazda fosternest. Alfabetycznie - jakie gniazda wypadają w najniższej kategorii wagowej? (Spróbuj wyświetlić tylko pierwsze wiersze posortowanego zestawu danych, reszty nie potrzebujesz - np. 10-20 pierwszych).
Spodziewany efekt
arrange(dane, weight, fosternest)[1:20,]## tarsus back dam fosternest hatchdate sex weight habitat bill_length
## 1 13.10 549.4057 R186912 B1102 54 Fem 8.6 park 10.467
## 2 12.57 551.1584 R186908 B1602 53 Fem 8.6 park 10.071
## 3 13.10 549.2661 R187539 A102 47 Fem 8.8 forest 10.416
## 4 13.25 548.5838 R186912 B1102 54 Fem 8.8 park 10.559
## 5 12.95 550.3171 R187523 E2002 46 Fem 8.8 park 10.332
## 6 13.10 551.0305 R046109 C402 44 Fem 8.9 park 10.448
## 7 13.25 550.9402 R187557 F1902 47 Male 8.9 forest 10.600
## 8 13.55 551.3757 R187557 F2102 47 Fem 8.9 park 10.885
## 9 13.40 549.0100 R187588 F2402 47 Fem 8.9 park 10.674
## 10 13.48 548.0596 K983388 H1302 48 Fem 8.9 forest 10.783
## 11 13.25 549.0567 Fem20 H1302 48 Fem 8.9 forest 10.558
## 12 13.40 549.0138 R187527 A1802 45 Fem 9.0 park 10.772
## 13 13.02 554.1662 R187562 A2302 47 Fem 9.0 forest 10.441
## 14 13.25 550.1252 R187562 A602 47 Fem 9.0 park 10.616
## 15 13.55 549.6995 R187902 B1102 54 UNK 9.0 forest 10.850
## 16 13.86 549.4298 R187595 C2402 48 Fem 9.0 forest 11.034
## 17 13.33 548.9120 R186907 E1002 53 Fem 9.0 park 10.648
## 18 13.86 548.7900 R187531 F1702 46 Fem 9.0 forest 11.106
## 19 13.48 550.5093 R187916 G1602 49 Fem 9.0 park 10.798
## 20 13.78 549.4879 R187931 G2202 49 Fem 9.0 park 11.008
## bill_depth
## 1 0.584
## 2 0.574
## 3 0.573
## 4 0.526
## 5 0.517
## 6 0.587
## 7 0.580
## 8 0.597
## 9 0.577
## 10 0.620
## 11 0.548
## 12 0.593
## 13 0.632
## 14 0.601
## 15 0.577
## 16 0.570
## 17 0.600
## 18 0.599
## 19 0.604
## 20 0.647
Zadanie Zmodyfikuj powyższe wyszukiwanie zmieniając
cechę na skok, oraz sortując malejąco (użyj modyfikatora
desc na zmiennej tarsus).
Spodziewany wynik
arrange(dane, desc(tarsus), fosternest)[1:20,]## tarsus back dam fosternest hatchdate sex weight habitat bill_length
## 1 16.81 548.3708 R187545 D1002 47 Male 11.4 park 13.399
## 2 15.98 549.7797 R187537 A2702 47 Male 10.8 park 12.828
## 3 15.67 550.5012 R187914 G602 49 Male 10.5 park 12.535
## 4 15.64 550.5361 R187292 E1802 47 Male 10.2 forest 12.471
## 5 15.60 550.4936 R187155 F2402 47 Male 11.2 park 12.477
## 6 15.52 551.3060 R186908 B1602 53 UNK 10.4 forest 12.391
## 7 15.52 549.4172 R186908 B1702 53 Male 10.4 park 12.418
## 8 15.52 550.8471 R187030 C1302 46 Male 10.5 forest 12.403
## 9 15.52 551.7863 R187590 D202 49 Male 10.5 forest 12.394
## 10 15.45 548.3911 R187030 C1302 46 Fem 10.2 forest 12.348
## 11 15.45 549.9689 R187030 C1302 46 Male 10.6 forest 12.324
## 12 15.45 550.7465 R187927 D402 50 Male 9.8 park 12.298
## 13 15.45 551.0630 R187577 E1802 47 Male 10.4 park 12.358
## 14 15.45 550.6196 R187521 E2002 46 Male 10.4 park 12.345
## 15 15.45 549.8546 R187399 E402 47 Male 10.4 park 12.350
## 16 15.45 552.1680 R187531 F1702 46 Male 10.4 forest 12.378
## 17 15.45 549.6147 R187030 F1702 46 Male 10.6 park 12.424
## 18 15.45 548.6120 R187030 F1702 46 Male 10.3 forest 12.384
## 19 15.45 547.8379 Fem20 H1302 48 Male 10.8 park 12.389
## 20 15.37 549.5868 R187515 A1802 45 Male 10.5 park 12.289
## bill_depth
## 1 0.797
## 2 0.690
## 3 0.711
## 4 0.709
## 5 0.696
## 6 0.693
## 7 0.674
## 8 0.697
## 9 0.725
## 10 0.675
## 11 0.706
## 12 0.671
## 13 0.682
## 14 0.624
## 15 0.669
## 16 0.669
## 17 0.675
## 18 0.696
## 19 0.728
## 20 0.633
\(\blacksquare\)
3. Selekcja zmiennych dzięki select()
select() pozwala wybierać pojedyncze kolumny z danych,
albo za pomocą ich pełnej nazwy albo za pomocą modyfikatorów
dopasowujących kolumny do zadanych warunków:
starts_with()- dopasowuje kolumna na podstawie początkowych znakówends_with()- podobnie, jak wyżej, ale dla znaków końcowychcontains()- dopasowuje na podstawie ciągu znaków zawartego gdzieś w nazwie kolumnymatches()- dopasowuje kolumny na podstawie wyrażenia regularnego (patrz linki na górze)
Wybierz z zestawu danych kolumny od 1 do 4 oraz 8 (we wszystkich
przykładach wyświetlamy pierwszych 20 wierszy za pomocą
head()):
head(select(dane, 1:4, 8),
20)## tarsus back dam fosternest habitat
## 1 13.55 551.3757 R187557 F2102 park
## 2 15.07 549.0884 R187559 F1902 forest
## 3 14.99 550.1739 R187568 A602 forest
## 4 14.69 550.3067 R187518 A1302 forest
## 5 14.46 549.6392 R187528 A2602 forest
## 6 13.93 551.8693 R187945 C2302 forest
## 7 13.93 549.4878 Fem3 C1902 park
## 8 15.45 548.3911 R187030 C1302 forest
## 9 14.31 550.0810 R187517 C602 park
## 10 14.46 549.8523 R187523 B2202 forest
## 11 14.99 547.7403 R186902 B1402 forest
## 12 14.08 548.9847 R187400 B1002 park
## 13 14.77 551.1920 R187932 B502 park
## 14 13.63 549.2743 R187582 D1202 forest
## 15 16.81 548.3708 R187545 D1002 park
## 16 14.24 551.1921 R187546 D902 park
## 17 14.92 551.5066 R187590 D202 forest
## 18 14.46 549.7494 R187548 E902 park
## 19 15.37 552.3874 R187594 E302 forest
## 20 13.40 549.0100 R187588 F2402 park
Wybierz kolumny zawierające pomiary strukturalne (skok oraz dwa pomiary dzioba):
head(select(dane, tarsus, starts_with("bill")),
20)## tarsus bill_length bill_depth
## 1 13.55 10.885 0.597
## 2 15.07 12.056 0.682
## 3 14.99 12.000 0.673
## 4 14.69 11.733 0.639
## 5 14.46 11.565 0.631
## 6 13.93 11.165 0.626
## 7 13.93 11.154 0.591
## 8 15.45 12.348 0.675
## 9 14.31 11.434 0.629
## 10 14.46 11.564 0.609
## 11 14.99 11.977 0.636
## 12 14.08 11.251 0.685
## 13 14.77 11.828 0.638
## 14 13.63 10.904 0.612
## 15 16.81 13.399 0.797
## 16 14.24 11.409 0.679
## 17 14.92 11.929 0.700
## 18 14.46 11.512 0.669
## 19 15.37 12.259 0.641
## 20 13.40 10.674 0.577
Przenieś kolumnę hatchdate na sam koniec tabeli:
head(select(dane, -hatchdate, everything(), hatchdate),
20)## tarsus back dam fosternest sex weight habitat bill_length
## 1 13.55 551.3757 R187557 F2102 Fem 8.9 park 10.885
## 2 15.07 549.0884 R187559 F1902 Male 10.5 forest 12.056
## 3 14.99 550.1739 R187568 A602 Male 9.9 forest 12.000
## 4 14.69 550.3067 R187518 A1302 Male 9.7 forest 11.733
## 5 14.46 549.6392 R187528 A2602 Fem 9.7 forest 11.565
## 6 13.93 551.8693 R187945 C2302 Fem 9.3 forest 11.165
## 7 13.93 549.4878 Fem3 C1902 Male 9.3 park 11.154
## 8 15.45 548.3911 R187030 C1302 Fem 10.2 forest 12.348
## 9 14.31 550.0810 R187517 C602 Fem 9.5 park 11.434
## 10 14.46 549.8523 R187523 B2202 Fem 9.6 forest 11.564
## 11 14.99 547.7403 R186902 B1402 Male 10.0 forest 11.977
## 12 14.08 548.9847 R187400 B1002 Fem 9.7 park 11.251
## 13 14.77 551.1920 R187932 B502 Male 9.9 park 11.828
## 14 13.63 549.2743 R187582 D1202 Male 9.2 forest 10.904
## 15 16.81 548.3708 R187545 D1002 Male 11.4 park 13.399
## 16 14.24 551.1921 R187546 D902 Fem 9.8 park 11.409
## 17 14.92 551.5066 R187590 D202 Fem 10.1 forest 11.929
## 18 14.46 549.7494 R187548 E902 Fem 9.6 park 11.512
## 19 15.37 552.3874 R187594 E302 Male 10.4 forest 12.259
## 20 13.40 549.0100 R187588 F2402 Fem 8.9 park 10.674
## bill_depth hatchdate
## 1 0.597 47
## 2 0.682 47
## 3 0.673 47
## 4 0.639 45
## 5 0.631 45
## 6 0.626 49
## 7 0.591 47
## 8 0.675 46
## 9 0.629 44
## 10 0.609 46
## 11 0.636 52
## 12 0.685 47
## 13 0.638 49
## 14 0.612 48
## 15 0.797 47
## 16 0.679 47
## 17 0.700 49
## 18 0.669 48
## 19 0.641 49
## 20 0.577 47
Zadanie Usuń z tabeli danych kolumny z pomiarami morfometrycznymi (skok, waga, dziób).
Spodziewany wynik
head(select(dane, -tarsus, weight, -starts_with("bill")),
20)## back dam fosternest hatchdate sex weight habitat
## 1 551.3757 R187557 F2102 47 Fem 8.9 park
## 2 549.0884 R187559 F1902 47 Male 10.5 forest
## 3 550.1739 R187568 A602 47 Male 9.9 forest
## 4 550.3067 R187518 A1302 45 Male 9.7 forest
## 5 549.6392 R187528 A2602 45 Fem 9.7 forest
## 6 551.8693 R187945 C2302 49 Fem 9.3 forest
## 7 549.4878 Fem3 C1902 47 Male 9.3 park
## 8 548.3911 R187030 C1302 46 Fem 10.2 forest
## 9 550.0810 R187517 C602 44 Fem 9.5 park
## 10 549.8523 R187523 B2202 46 Fem 9.6 forest
## 11 547.7403 R186902 B1402 52 Male 10.0 forest
## 12 548.9847 R187400 B1002 47 Fem 9.7 park
## 13 551.1920 R187932 B502 49 Male 9.9 park
## 14 549.2743 R187582 D1202 48 Male 9.2 forest
## 15 548.3708 R187545 D1002 47 Male 11.4 park
## 16 551.1921 R187546 D902 47 Fem 9.8 park
## 17 551.5066 R187590 D202 49 Fem 10.1 forest
## 18 549.7494 R187548 E902 48 Fem 9.6 park
## 19 552.3874 R187594 E302 49 Male 10.4 forest
## 20 549.0100 R187588 F2402 47 Fem 8.9 park
Zadanie Pozostaw w tabeli tylko kolumny zawierające w nazwie przynajmniej jedno “e”.
Spodziewany wynik
head(select(dane, contains("e")),
20)## fosternest hatchdate sex weight bill_length bill_depth
## 1 F2102 47 Fem 8.9 10.885 0.597
## 2 F1902 47 Male 10.5 12.056 0.682
## 3 A602 47 Male 9.9 12.000 0.673
## 4 A1302 45 Male 9.7 11.733 0.639
## 5 A2602 45 Fem 9.7 11.565 0.631
## 6 C2302 49 Fem 9.3 11.165 0.626
## 7 C1902 47 Male 9.3 11.154 0.591
## 8 C1302 46 Fem 10.2 12.348 0.675
## 9 C602 44 Fem 9.5 11.434 0.629
## 10 B2202 46 Fem 9.6 11.564 0.609
## 11 B1402 52 Male 10.0 11.977 0.636
## 12 B1002 47 Fem 9.7 11.251 0.685
## 13 B502 49 Male 9.9 11.828 0.638
## 14 D1202 48 Male 9.2 10.904 0.612
## 15 D1002 47 Male 11.4 13.399 0.797
## 16 D902 47 Fem 9.8 11.409 0.679
## 17 D202 49 Fem 10.1 11.929 0.700
## 18 E902 48 Fem 9.6 11.512 0.669
## 19 E302 49 Male 10.4 12.259 0.641
## 20 F2402 47 Fem 8.9 10.674 0.577
\(\blacksquare\)
4. Tworzenie nowych zmiennych dzięki mutate()
Stwórz zmienną bill_ratio będącą stosunkiem długości dzioba
do jego głębokości. Użyj do tego skróconego zestawu dane2,
w którym pozostawisz tylko kolumny morfometryczne oraz ID matki
(dam) oraz gniazdo (fosternest):
dane2 <- select(dane, dam, fosternest, tarsus, weight, starts_with("bill"))
head(mutate(dane2, bill_ratio = round(bill_length/bill_depth, 3)),
20)## dam fosternest tarsus weight bill_length bill_depth bill_ratio
## 1 R187557 F2102 13.55 8.9 10.885 0.597 18.233
## 2 R187559 F1902 15.07 10.5 12.056 0.682 17.677
## 3 R187568 A602 14.99 9.9 12.000 0.673 17.831
## 4 R187518 A1302 14.69 9.7 11.733 0.639 18.362
## 5 R187528 A2602 14.46 9.7 11.565 0.631 18.328
## 6 R187945 C2302 13.93 9.3 11.165 0.626 17.835
## 7 Fem3 C1902 13.93 9.3 11.154 0.591 18.873
## 8 R187030 C1302 15.45 10.2 12.348 0.675 18.293
## 9 R187517 C602 14.31 9.5 11.434 0.629 18.178
## 10 R187523 B2202 14.46 9.6 11.564 0.609 18.989
## 11 R186902 B1402 14.99 10.0 11.977 0.636 18.832
## 12 R187400 B1002 14.08 9.7 11.251 0.685 16.425
## 13 R187932 B502 14.77 9.9 11.828 0.638 18.539
## 14 R187582 D1202 13.63 9.2 10.904 0.612 17.817
## 15 R187545 D1002 16.81 11.4 13.399 0.797 16.812
## 16 R187546 D902 14.24 9.8 11.409 0.679 16.803
## 17 R187590 D202 14.92 10.1 11.929 0.700 17.041
## 18 R187548 E902 14.46 9.6 11.512 0.669 17.208
## 19 R187594 E302 15.37 10.4 12.259 0.641 19.125
## 20 R187588 F2402 13.40 8.9 10.674 0.577 18.499
Zadanie Za pomocą funkcji paste()
stwórz zmienną zawierającą kombinację ID matki (dam) i
gniazda wychowywania (fosternest). Nazwij tą zmienną
“crossfoster”.
Spodziewany wynik
head(mutate(dane2, crossfoster = paste(dam, fosternest, sep = "_")),
20)## dam fosternest tarsus weight bill_length bill_depth crossfoster
## 1 R187557 F2102 13.55 8.9 10.885 0.597 R187557_F2102
## 2 R187559 F1902 15.07 10.5 12.056 0.682 R187559_F1902
## 3 R187568 A602 14.99 9.9 12.000 0.673 R187568_A602
## 4 R187518 A1302 14.69 9.7 11.733 0.639 R187518_A1302
## 5 R187528 A2602 14.46 9.7 11.565 0.631 R187528_A2602
## 6 R187945 C2302 13.93 9.3 11.165 0.626 R187945_C2302
## 7 Fem3 C1902 13.93 9.3 11.154 0.591 Fem3_C1902
## 8 R187030 C1302 15.45 10.2 12.348 0.675 R187030_C1302
## 9 R187517 C602 14.31 9.5 11.434 0.629 R187517_C602
## 10 R187523 B2202 14.46 9.6 11.564 0.609 R187523_B2202
## 11 R186902 B1402 14.99 10.0 11.977 0.636 R186902_B1402
## 12 R187400 B1002 14.08 9.7 11.251 0.685 R187400_B1002
## 13 R187932 B502 14.77 9.9 11.828 0.638 R187932_B502
## 14 R187582 D1202 13.63 9.2 10.904 0.612 R187582_D1202
## 15 R187545 D1002 16.81 11.4 13.399 0.797 R187545_D1002
## 16 R187546 D902 14.24 9.8 11.409 0.679 R187546_D902
## 17 R187590 D202 14.92 10.1 11.929 0.700 R187590_D202
## 18 R187548 E902 14.46 9.6 11.512 0.669 R187548_E902
## 19 R187594 E302 15.37 10.4 12.259 0.641 R187594_E302
## 20 R187588 F2402 13.40 8.9 10.674 0.577 R187588_F2402
\(\blacksquare\)
5. Grupowanie obserwacji przez group_by()
Grupowanie pozwala stworzyć w zestawie danych “podzbiory” traktowane jako oddzielne pule przy liczeniu określonych podsumowań.
Porównaj taką operację:
head(mutate(dane2, mean_tars = mean(tarsus)),
20)## dam fosternest tarsus weight bill_length bill_depth mean_tars
## 1 R187557 F2102 13.55 8.9 10.885 0.597 14.5005
## 2 R187559 F1902 15.07 10.5 12.056 0.682 14.5005
## 3 R187568 A602 14.99 9.9 12.000 0.673 14.5005
## 4 R187518 A1302 14.69 9.7 11.733 0.639 14.5005
## 5 R187528 A2602 14.46 9.7 11.565 0.631 14.5005
## 6 R187945 C2302 13.93 9.3 11.165 0.626 14.5005
## 7 Fem3 C1902 13.93 9.3 11.154 0.591 14.5005
## 8 R187030 C1302 15.45 10.2 12.348 0.675 14.5005
## 9 R187517 C602 14.31 9.5 11.434 0.629 14.5005
## 10 R187523 B2202 14.46 9.6 11.564 0.609 14.5005
## 11 R186902 B1402 14.99 10.0 11.977 0.636 14.5005
## 12 R187400 B1002 14.08 9.7 11.251 0.685 14.5005
## 13 R187932 B502 14.77 9.9 11.828 0.638 14.5005
## 14 R187582 D1202 13.63 9.2 10.904 0.612 14.5005
## 15 R187545 D1002 16.81 11.4 13.399 0.797 14.5005
## 16 R187546 D902 14.24 9.8 11.409 0.679 14.5005
## 17 R187590 D202 14.92 10.1 11.929 0.700 14.5005
## 18 R187548 E902 14.46 9.6 11.512 0.669 14.5005
## 19 R187594 E302 15.37 10.4 12.259 0.641 14.5005
## 20 R187588 F2402 13.40 8.9 10.674 0.577 14.5005
z taką:
dane3 <- group_by(dane2, dam)
head(mutate(dane3, mean_tars = mean(tarsus)),
20)## # A tibble: 20 × 7
## # Groups: dam [20]
## dam fosternest tarsus weight bill_length bill_depth mean_tars
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 R187557 F2102 13.6 8.9 10.9 0.597 13.9
## 2 R187559 F1902 15.1 10.5 12.1 0.682 14.8
## 3 R187568 A602 15.0 9.9 12 0.673 14.3
## 4 R187518 A1302 14.7 9.7 11.7 0.639 14.5
## 5 R187528 A2602 14.5 9.7 11.6 0.631 14.4
## 6 R187945 C2302 13.9 9.3 11.2 0.626 14.7
## 7 Fem3 C1902 13.9 9.3 11.2 0.591 14.5
## 8 R187030 C1302 15.4 10.2 12.3 0.675 15.2
## 9 R187517 C602 14.3 9.5 11.4 0.629 14.8
## 10 R187523 B2202 14.5 9.6 11.6 0.609 14.2
## 11 R186902 B1402 15.0 10 12.0 0.636 14.2
## 12 R187400 B1002 14.1 9.7 11.3 0.685 14.7
## 13 R187932 B502 14.8 9.9 11.8 0.638 14.2
## 14 R187582 D1202 13.6 9.2 10.9 0.612 14.5
## 15 R187545 D1002 16.8 11.4 13.4 0.797 14.7
## 16 R187546 D902 14.2 9.8 11.4 0.679 14.6
## 17 R187590 D202 14.9 10.1 11.9 0.7 15.0
## 18 R187548 E902 14.5 9.6 11.5 0.669 14.2
## 19 R187594 E302 15.4 10.4 12.3 0.641 14.7
## 20 R187588 F2402 13.4 8.9 10.7 0.577 14.3
Aby lepiej zobaczyć co właściwie się tutaj stało - posortujmy wynikową tabelę wg ID matki:
dane3 <- group_by(dane2, dam)
dane3 <- mutate(dane3, mean_tars = mean(tarsus))
arrange(dane3, dam)## # A tibble: 828 × 7
## # Groups: dam [106]
## dam fosternest tarsus weight bill_length bill_depth mean_tars
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Fem2 H702 14.5 10 11.6 0.645 14.4
## 2 Fem2 H702 13.8 9.5 11.1 0.619 14.4
## 3 Fem2 A18B02 14.4 9.5 11.5 0.634 14.4
## 4 Fem2 A18B02 14.7 9.5 11.8 0.658 14.4
## 5 Fem2 A18B02 14.8 9.8 11.8 0.607 14.4
## 6 Fem20 H1102 13.9 9.4 11.1 0.567 14.5
## 7 Fem20 H1102 14.5 9.9 11.6 0.665 14.5
## 8 Fem20 H1102 14.8 9.9 11.9 0.643 14.5
## 9 Fem20 H1102 14.1 9.4 11.2 0.591 14.5
## 10 Fem20 H1102 14.6 9.8 11.7 0.643 14.5
## # … with 818 more rows
Czy zauważyłeś zmianę tabeli danych na nowy typ obiektu zwany
tibble?
Wykonanie powyższej serii operacji bardzo ułatwia użycie tzw. operatora pipe, pozwalającego łączyć dłuższe ciągi operacji w jeden spójny “potok”:
group_by(dane2, dam) %>%
mutate(mean_tars = mean(tarsus)) %>%
arrange(dam)## # A tibble: 828 × 7
## # Groups: dam [106]
## dam fosternest tarsus weight bill_length bill_depth mean_tars
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Fem2 H702 14.5 10 11.6 0.645 14.4
## 2 Fem2 H702 13.8 9.5 11.1 0.619 14.4
## 3 Fem2 A18B02 14.4 9.5 11.5 0.634 14.4
## 4 Fem2 A18B02 14.7 9.5 11.8 0.658 14.4
## 5 Fem2 A18B02 14.8 9.8 11.8 0.607 14.4
## 6 Fem20 H1102 13.9 9.4 11.1 0.567 14.5
## 7 Fem20 H1102 14.5 9.9 11.6 0.665 14.5
## 8 Fem20 H1102 14.8 9.9 11.9 0.643 14.5
## 9 Fem20 H1102 14.1 9.4 11.2 0.591 14.5
## 10 Fem20 H1102 14.6 9.8 11.7 0.643 14.5
## # … with 818 more rows
Zadanie Korzystając z operatora %>%
zmodyfikuj powyższy ciąg tak, by na końcu powstała tabela o
następujących właściwościach:
- grupowanie ze względu na
fosternest; - obliczenie średniej wagi piskląt w obrębie każdego gniazda wychowywania;
- usunięcie z tabeli zmiennych niepotrzebnych (cechy dzioba)
- obliczenie stosunku wagi średniej w danym gnieździe do średniej wagi
w całej populacji. (Uwaga! Ten etap wymaga usunięcia grupowania
za pomocą funkcji
ungroup())
Aby wyświetlić więcej niż 10 wierszy tabeli dodaj do potoku funkcji
print() z opcją n = 30 (jeśli chcesz
wyświetlić 30 wierszy). Dlaczego w trakcie obliczeń musimy zastosować
ungroup()? Co się stanie, jeśli tego nie zrobimy?
Spodziewany wynik
group_by(dane2, fosternest) %>%
mutate(mean_weight = mean(weight)) %>%
select(-starts_with("bill")) %>%
arrange(fosternest) %>%
ungroup() %>%
mutate(mean_ratio = mean_weight/mean(weight)) %>%
print(n = 30)## # A tibble: 828 × 6
## dam fosternest tarsus weight mean_weight mean_ratio
## <fct> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 R187537 A1002 14.7 10 9.64 0.985
## 2 R187537 A1002 14.8 9.8 9.64 0.985
## 3 R187537 A1002 14.9 9.7 9.64 0.985
## 4 R187537 A1002 14.5 9.6 9.64 0.985
## 5 R187563 A1002 13.2 9.3 9.64 0.985
## 6 R187563 A1002 13.6 9.2 9.64 0.985
## 7 R187563 A1002 14.0 9.2 9.64 0.985
## 8 R187563 A1002 14.7 10.1 9.64 0.985
## 9 R187563 A1002 14.4 9.9 9.64 0.985
## 10 R187539 A102 14.2 9.8 9.67 0.987
## 11 R187539 A102 14.1 9.5 9.67 0.987
## 12 R187539 A102 15.1 10.1 9.67 0.987
## 13 R187539 A102 13.1 8.8 9.67 0.987
## 14 R187566 A102 14.6 9.4 9.67 0.987
## 15 R187566 A102 14.7 9.9 9.67 0.987
## 16 R187566 A102 14.4 9.6 9.67 0.987
## 17 R187566 A102 14.5 9.8 9.67 0.987
## 18 R187566 A102 15.1 10.1 9.67 0.987
## 19 R186910 A1202 14.5 10.2 9.7 0.991
## 20 R186910 A1202 13.5 9.1 9.7 0.991
## 21 R186910 A1202 14.8 10.4 9.7 0.991
## 22 R186910 A1202 14.1 9.2 9.7 0.991
## 23 R186910 A1202 15.1 10.3 9.7 0.991
## 24 R186910 A1202 14.1 9.5 9.7 0.991
## 25 R186910 A1202 13.9 9.2 9.7 0.991
## 26 R187518 A1302 14.7 9.7 9.84 1.00
## 27 R187518 A1302 14.8 10.1 9.84 1.00
## 28 R187518 A1302 14.1 9.4 9.84 1.00
## 29 R187518 A1302 14.6 9.9 9.84 1.00
## 30 R187518 A1302 14.1 9.7 9.84 1.00
## # … with 798 more rows
\(\blacksquare\)
6. Podsumowanie z summarise()
Przefiltruj dane tak, by zostały w nich tylko ptaki ze
znaną płcią (sex == "Male" lub sex == "Fem") i
podsumuj je wyliczając dla obserwacji średnią długość skoku, liczbę
obserwacji oraz błąd standardowy skoku (liczony jako
SD/sqrt(N)):
dane %>%
filter(sex != "UNK") %>%
summarise(mean_t = mean(tarsus),
N = n(),
se_t = sd(tarsus)/sqrt(n()))## mean_t N se_t
## 1 14.50553 781 0.01797891
Zadanie Powtórz powyższe obliczenia, ale dla ptaków
w każdej możliwej kombinacji sex oraz habitat
(tzn. samce w lesie, samce w parku, etc.)
Spodziewany wynik
dane %>%
filter(sex != "UNK") %>%
group_by(sex, habitat) %>%
summarise(mean_t = mean(tarsus),
N = n(),
se_t = sd(tarsus)/sqrt(n()))## `summarise()` has grouped output by 'sex'. You can override using the `.groups`
## argument.
## # A tibble: 4 × 5
## # Groups: sex [2]
## sex habitat mean_t N se_t
## <fct> <fct> <dbl> <int> <dbl>
## 1 Fem forest 14.3 188 0.0334
## 2 Fem park 14.3 185 0.0343
## 3 Male forest 14.7 205 0.0325
## 4 Male park 14.7 203 0.0325
\(\blacksquare\)
7. Szerokie i długie zestawy danych
Oryginalny zestaw dane jest typu “długiego”: każda
obserwacja to jeden osobnik, a zmienne kategoryczne zajmują po jednej
kolumnie zawierającej przypisania każdego osobnika do określonej
kategorii. Często jednak dane, które dostajemy są typu “szerokiego” (np.
grupy eksperymentalne czy płcie znajdują się w osobnych kolumnach).
Między tymi dwoma typami danych można sprawnie się przemieszczać za
pomocą funkcji pivot_wider() i
pivot_longer().
By dowiedzieć się o nich nieco więcej - zajrzyj tutaj:
vignette("pivot")Stwórzmy przykładowy zestaw danych zawierający średnie ciężary
piskląt i ich wariancje według płci, pogrupowane według
fosternest:
dane_l <- dane %>%
filter(sex != "UNK") %>%
group_by(fosternest, sex) %>%
summarise(mean_w = mean(weight), var_w = var(weight))## `summarise()` has grouped output by 'fosternest'. You can override using the
## `.groups` argument.
Taki zestaw danych w formie długiej możemy przerobić na szeroki tworząc osobne kolumny dla samców i samic w danym gnieździe.
dane_l %>%
pivot_wider(names_from = sex, values_from = mean_w)## # A tibble: 201 × 4
## # Groups: fosternest [104]
## fosternest var_w Fem Male
## <fct> <dbl> <dbl> <dbl>
## 1 A1002 0.139 9.57 NA
## 2 A1002 0.07 NA 9.8
## 3 A102 0.249 9.48 NA
## 4 A102 0.0770 NA 9.82
## 5 A1202 0.0433 9.27 NA
## 6 A1202 0.00500 NA 10.4
## 7 A1302 0 9.7 NA
## 8 A1302 0.0817 NA 9.88
## 9 A1602 0.0960 9.5 NA
## 10 A1602 0.180 NA 9.7
## # … with 191 more rows
Powyższa komenda nie zadziałała tak, jak chcemy, bo pozostawiliśmy w
danych kolumnę z wariancją, zawierającą unikatowe wpisy “dublujące”
wiersze z każdym fosternest. Musimy albo się tej kolumny
pozbyć, albo też uwzględnić ją w pivot_wider:
dane_l %>%
select(-var_w) %>%
pivot_wider(names_from = sex, values_from = mean_w)## # A tibble: 104 × 3
## # Groups: fosternest [104]
## fosternest Fem Male
## <fct> <dbl> <dbl>
## 1 A1002 9.57 9.8
## 2 A102 9.48 9.82
## 3 A1202 9.27 10.4
## 4 A1302 9.7 9.88
## 5 A1602 9.5 9.7
## 6 A1802 9.6 9.95
## 7 A18B02 10 9.72
## 8 A2202 9.48 10
## 9 A22B02 9.33 9.48
## 10 A2302 9.57 9.7
## # … with 94 more rows
dane_l %>%
pivot_wider(names_from = sex, values_from = c(mean_w, var_w))## # A tibble: 104 × 5
## # Groups: fosternest [104]
## fosternest mean_w_Fem mean_w_Male var_w_Fem var_w_Male
## <fct> <dbl> <dbl> <dbl> <dbl>
## 1 A1002 9.57 9.8 0.139 0.07
## 2 A102 9.48 9.82 0.249 0.0770
## 3 A1202 9.27 10.4 0.0433 0.00500
## 4 A1302 9.7 9.88 0 0.0817
## 5 A1602 9.5 9.7 0.0960 0.180
## 6 A1802 9.6 9.95 0.260 0.605
## 7 A18B02 10 9.72 0.0200 0.0870
## 8 A2202 9.48 10 0.192 0.167
## 9 A22B02 9.33 9.48 0.0233 0.0692
## 10 A2302 9.57 9.7 0.171 0.07
## # … with 94 more rows
Zapiszmy dane w nowym obiekcie dane_w.
dane_w <- dane_l %>%
select(-var_w) %>%
pivot_wider(names_from = sex, values_from = mean_w)Zadanie Spróbuj odwrócić działanie zamieniając
“szeroki” styl danych na “długi” - korzystając z funkcji
pivot_longer(). Nazwij zmienną “płeć”
animal_sex, a zmienną z ciężarem ciała - ’body_weight`.
Spodziewany wynik
dane_w %>%
pivot_longer(cols = c(Fem, Male),
names_to = "animal_sex",
values_to = "body_weight")## # A tibble: 208 × 3
## # Groups: fosternest [104]
## fosternest animal_sex body_weight
## <fct> <chr> <dbl>
## 1 A1002 Fem 9.57
## 2 A1002 Male 9.8
## 3 A102 Fem 9.48
## 4 A102 Male 9.82
## 5 A1202 Fem 9.27
## 6 A1202 Male 10.4
## 7 A1302 Fem 9.7
## 8 A1302 Male 9.88
## 9 A1602 Fem 9.5
## 10 A1602 Male 9.7
## # … with 198 more rows
B. Prosta statystyka
1. Czy samce są większe od samic?
W R proste testy statystyczne wykonujemy za pomocą serii
funkcji z rodziny __.test - np. takich:
| Funkcja | Test | Zastosowanie |
|---|---|---|
t.test() |
test t-Studenta | porównywanie średnich |
var.test() |
test F | porównywanie wariancji (zmienności) |
cor.test() |
test korelacji | testowanie siły związku między zmiennymi |
chisq.test() |
test \(\chi^2\) | testowanie zmiennych kategorycznych i tabel wielopolowych |
kruskal.test() |
test Kruskall-Wallisa | nieparametryczne porównywanie grup |
wilcox.test |
test Wilcoxona | nieparametryczny test par wiązanych |
Sprawdźmy za pomocą testu korelacji czy waga piskląt związana jest istotnie z ich długością skoku. Jaka jest: wartość p testu, obliczona korelacja i jej 95% przedział ufności? Czy odrzucasz hipotezę zerową? Jaka ona jest?
test1 <- cor.test(dane$tarsus, dane$weight)
test1##
## Pearson's product-moment correlation
##
## data: dane$tarsus and dane$weight
## t = 48.044, df = 826, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8391000 0.8751343
## sample estimates:
## cor
## 0.8581706
test1$p.value## [1] 2.096071e-241
Prosty wykres rozrzutu dobrze obrazuje otrzymany wynik:
plot(dane$tarsus, dane$weight,
xlab = "Weight", ylab = "Tarsus")Zadanie Wiedząc, że opcja wykresu col
może m.in. przyjmować wartości będące zmiennymi kategorycznymi - jak
spróbujesz dodać do wykresy informację o płci, zmieniającą kolor
obserwacji w zależności od płci? Wskazówka Zamiast
dane$zmienna - użyj
jitter(dane$zmienna, 3).
Sugerowany wynik
plot(jitter(dane$tarsus, 3), jitter(dane$weight, 3),
xlab = "Weight", ylab = "Tarsus",
col = dane$sex)Zadanie Z wykresu jasno wynika, że samce powinny być
większe od samic pod względem zarówno masy, jak i skoku. Przetestujmy tą
różnicę za pomocą testu t-Studenta. Wykonaj go ekstrahując z
danych skok samic, a potem samców (za pomocą funkcji
subset(), lub w inny znany Ci sposób). Dla przypomnienia:
subset(dane, habitat == "forest")$weight pozwoliłby
wyciągnąć z danych wagi ptaków żyjących w lesie. Wykonaj go w podobny
sposób, jak cor.test(). Jaka jest Twoja decyzja: odrzucasz
czy przyjmujesz hipotezę zerową o braku różnic między płciami?
Sugerowany wynik
test2 <- t.test(subset(dane, sex == "Male")$tarsus,
subset(dane, sex == "Fem")$tarsus)
test2##
## Welch Two Sample t-test
##
## data: subset(dane, sex == "Male")$tarsus and subset(dane, sex == "Fem")$tarsus
## t = 11.759, df = 773.68, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.3250593 0.4553372
## sample estimates:
## mean of x mean of y
## 14.69189 14.30169
W przypadku danych w formacie “długim” posługiwanie się funkcją
t.test() jest znacznie łatwiejsze w konwencji tzw. formuły.
Formuła to rodzaj wyrażenia matematycznego, które pozwala
zapisać zależności istniejące między zmienną/zmiennymi zależną/zależnymi
i niezależną/niezależnymi. W przypadku prostej zależności, gdzie mamy
jedną zmienną zależną y (np. skok) i jedną zmienną
niezależną x (np. płeć) - formułę możemy zapisać jako
y ~ x.
Zadanie Powtórz powyższy test t-Studenta
używają notacji formuły zamiast dwóch osobnych wektorów zawierających
dane samic i samców. Aby konwencja ta zadziałała będziesz potrzebować
wyrażenia opisującego zależność skoku od płci
(tarsus ~ sex) oraz dodatkowej opcji data w
t.test() wskazującej na to z jakiego obiektu pochodzą dane
(zajrzyj do ?t.test() w razie wątpliwości).
Uwaga! W dalszym ciągu musisz usunąć z danych
obserwacje sex == "UNK" - w teście t-Studenta
możemy porównywać tylko 2 grupy.
Sugerowany wynik
test3 <- t.test(tarsus ~ sex, data = subset(dane, sex != "UNK"))
test3##
## Welch Two Sample t-test
##
## data: tarsus by sex
## t = -11.759, df = 773.68, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Fem and group Male is not equal to 0
## 95 percent confidence interval:
## -0.4553372 -0.3250593
## sample estimates:
## mean in group Fem mean in group Male
## 14.30169 14.69189
Zadanie Przeprowadzany test t-Studenta
zakłada, że wariancje porównywanych grupo są różne. Przeprowadzając test
funkcją var.test() zweryfikuj to założenie, testując
hipotezę zerową, że wariancje są równe. Składnia testu jest bardzo
podobna do tej użytej powyżej z wykorzystaniem formuły.
Sugerowany wynik
test4 <- var.test(tarsus ~ sex, data = subset(dane, sex != "UNK"))
test4##
## F test to compare two variances
##
## data: tarsus by sex
## F = 0.98626, num df = 372, denom df = 407, p-value = 0.8928
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.808562 1.204097
## sample estimates:
## ratio of variances
## 0.9862598
\(\blacksquare\)
2. Czy samic lub samców jest więcej w lesie/w parku?
Aby odpowiedzieć na to pytanie - potrzebujemy podsumowania
zliczającego samce i samice w dwóch rodzajach siedliska. Najłatwiej
uzyskać takie podsumowanie stosując funkcję table():
dane2 <- filter(dane, sex != "UNK") %>% gdata::drop.levels()
sex_hab <- table(dane2$sex, dane2$habitat)
sex_hab##
## forest park
## Fem 188 185
## Male 205 203
Zadanie Wiedząc, że funkcja chisq.test
potrafi jako jeden z argumentów przyjmować całą tabelę wielopolową -
przeprowadź test chi-kwadrat hipotezy badającej zależność występowania
ptaków w 2 siedliskach w zależności od płci. Czy zależność ta jest
istotna? Co znajduje się w zbiorze stdres powstałego
obiektu?
Spodziewany wynik
test5 <- chisq.test(sex_hab)
test5##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: sex_hab
## X-squared = 0, df = 1, p-value = 1
test5$stdres##
## forest park
## Fem 0.04384568 -0.04384568
## Male -0.04384568 0.04384568
\(\blacksquare\)
3. Czy długość skoku zależy od daty wyklucia? Regresja prosta
Sprawdźmy to najpierw wizualnie:
plot(tarsus ~ hatchdate, data = dane)Nie wygląda to na zbyt silną zależność - możemy jednak sprawdzić to formalnie za pomocą modelu liniowego w najprostszej jego wersji. Model liniowy, jaki będziemy testować wygląda następująco:
\[y_i = a + bx_i + e_i\] czyli
\[tarsus_i = a + b*hatchdate_i +
e_i\] Do jego obliczenia służy funkcja lm()
(linear model) - co zaskakujące
przyjmuje ona argumenty w sposób bardzo podobny do np.
t.test() (korzystając z konwencji formuły).
Zadanie Stosując analogię wywołania funkcji
t.test(), użyj lm() by przetestować zależność
długości skoku sikory modrej od daty klucia
(tarsus ~ hatchdate). W pierwszym kroku stwórz obiekt
model zawierający wynik działania funkcji
lm(). W drugim kroku użyj summary() by
podejrzeć wynik działania obliczeń.
Sugerowany wynik
model <- lm(tarsus ~ hatchdate, data = dane)
summary(model)##
## Call:
## lm(formula = tarsus ~ hatchdate, data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8746 -0.3471 0.0468 0.3011 2.2925
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.08869 0.33443 45.117 <2e-16 ***
## hatchdate -0.01215 0.00690 -1.761 0.0786 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4996 on 826 degrees of freedom
## Multiple R-squared: 0.003741, Adjusted R-squared: 0.002535
## F-statistic: 3.102 on 1 and 826 DF, p-value: 0.07858
Przyjrzyj się wykresom, które powstają dzięki wywołaniu funkcji
plot(model)Pozwalają one ocenić, czy nasz model statystyczny spełnia założenia regresji prostej. Który z wykresów najlepiej oddaje test wizualny założenia zakładającego normalność wartości resztkowych modelu?
\(\blacksquare\)
C. Expressem przez bootstrapping
Do tej pory do naszych analiz wykorzystywaliśmy metody parametryczne
- możemy jednak spróbować wykorzystać inną, prostszą metodę, która
pozwoli nam uwolnić się od ograniczeń rozkładowych naszych danych.
Spójrzmy na przykład na rozkład zmiennej hatchdate - nie
jest on najpiękniejszy i na pierwszy rzut oka nie przypomina on
zgrabnego rozkładu normalnego:
hist(dane$hatchdate, 15)Spróbujmy skonstruować test oparty na bootstrappingu, którym
sprawdzimy czy ptaki z lasu (forest) różnią się pod kątem
daty klucia od ptaków z parków (park). Szybki rzut oka na
wykres pudełkowy nie wskazuje na zbyt wielkie różnice.
plot(hatchdate ~ habitat, data = dane)- Do skonstruowania bootstrappingu wykorzystamy pętlę
for, która pozwala powtarzać jakąś operację wiele razy. Jeśli chcielibyśmy np. wydrukować na ekranie (print()) jakiś tekst pięć razy - moglibyśmy osiągnąć to w następujący sposób:
for(i in 1:5) {
print("Tekst najbardziej oryginalny!")
}## [1] "Tekst najbardziej oryginalny!"
## [1] "Tekst najbardziej oryginalny!"
## [1] "Tekst najbardziej oryginalny!"
## [1] "Tekst najbardziej oryginalny!"
## [1] "Tekst najbardziej oryginalny!"
Zauważ, że pętla for wykorzystuje jakąś zmienną (tutaj
i), zmieniając jej wartość w zadanym zakresie (tutaj
1:5) i za każdym razem (dla każ∂ego i)
wykonując jakieś działania ujęte wewnątrz {}.
- Potrzebna nam będzie również metoda losowego wybierania pewnych
wierszy z naszych oryginalnych danych - ze zwracaniem. Idealnie do tego
nadaje się funkcja
sample(). Zobaczmy, jak losuje ona zadaną liczbę wartości (10) z podanej listy (liczby 1 do 20), ze zwracaniem (replace = T):
sample(1:20, size = 10, replace = T)## [1] 13 1 2 15 10 16 19 13 1 15
Uwaga! Funkcja ta używa generatora liczb pseudolosowych - więc u każdego wynik tej operacji będzie inny!
- Wreszcie - musimy mieć naszą statystykę testową. Niech będzie nią różnica między średnią datą klucia w sikor leśnych i parkowych \(d_0\) - jeśli wynosi ona zero, ptaki mają identyczne daty klucia:
d0 <- mean(filter(dane, habitat == "forest")$hatchdate) -
mean(filter(dane, habitat == "park")$hatchdate)
d0## [1] -0.1695764
Ujemna różnica wskazuje na nieco wyższą wartość (późniejsze klucie) w parku - pytanie, czy istotnie późniejsze?
Oto co powinna robić nasza procedura bootstrappingu:
- Krok 1: wylosuj z zestawu danych k = 828 obserwacji (wierszy) ze zwracaniem (tyle wierszy składa się na oryginalny zestaw danych).
- Krok 2: dla każdego z tych wylosowanych zestawów policz jego własną
statystykę
d- różnicę w średniej dacie klucia między lasem i parkiem. - Krok 3: Powtórz całą procedurę N = 1000 razy (jest to arbitralna
wielkość, im więcej - tym lepiej). Za każdym razem zapisz wynik do
zbiorczego obiekty
out. - Krok 4: Skonstruuj rozkład statystyki
di zobacz gdzie wypada w niej jej oryginalna (oparta na wejściowych danych) wartośćd0.
N <- 1000 # liczba zakładanych powtórzeń
k <- nrow(dane) # wielkość oryginalnego zbioru danych
# oryginalna statystyka testowa
d0 <- mean(filter(dane, habitat == "forest")$hatchdate) -
mean(filter(dane, habitat == "park")$hatchdate)
# obiekt wynikowy przygotowany na N=1000 końcowych wartości
out <- numeric(N)
for(i in 1:N) {
# losujemy 828 wierszy, ze zwracaniem
# a oryginalnej puli unikatowych 828 wierszy danych
rows <- sample(1:k, size = k, replace = T)
# tworzymy tymczasowe dane wewnątrz pętli
# zawierające tylko wylosowane wiersze
temp_dane <- dane[rows, ]
# oblicza na podstawie tak przygotowanych danych
# tymczasową statystykę testową dla powtórzenia
# nr i pętli
d_temp <- mean(filter(temp_dane, habitat == "forest")$hatchdate) -
mean(filter(temp_dane, habitat == "park")$hatchdate)
# zapisujemy wyliczoną statystykę
# do zmiennej out w miejscu i
out[i] <- d_temp
}Wyświetlmy rozkład symulowanej statystyki d oraz
położenie d0 na rozkładzie:
hist(out, 50)
abline(v = d0, lwd = 2)
abline(v = 0, col = "red", lwd = 2, lty = 2)Jak podejdziesz to odpowiedzi na pytanie - czy ta symulacja wskazuje na istotną statystycznie różnicę między symulowaną statystyką a wartością zero?